Guessing, if there is a Heart Attack or not
Statistical Methods and Laboratory, Part II - The Final Project
Sapienza Jeremy 1960498
6/18/2021
\[~\]
\[~\]
Introduction
\[~\]
Nowadays, with the improvement of the technologies there are an amount of data that allows us to understand if there is any problem or not in the healthy of a person.
Here in this project we explain and we try to predict whether a patient should be diagnosed with Heart Disease or not.
\[~\]
The Dataset
\[~\]
\[~\]
Kaggle, is the main platform where I found this interesting dataset where applying our main bayesian inference as the main scope of this project.
This dataset is composed by these features:
##
## -- Column specification --------------------------------------------------------
## cols(
## age = col_double(),
## sex = col_double(),
## cp = col_double(),
## trtbps = col_double(),
## chol = col_double(),
## fbs = col_double(),
## restecg = col_double(),
## thalachh = col_double(),
## exng = col_double(),
## oldpeak = col_double(),
## slp = col_double(),
## caa = col_double(),
## thall = col_double(),
## output = col_double()
## )## # A tibble: 6 x 14
## age sex cp trtbps chol fbs restecg thalachh exng oldpeak slp
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 63 1 3 145 233 1 0 150 0 2.3 0
## 2 37 1 2 130 250 0 1 187 0 3.5 0
## 3 41 0 1 130 204 0 0 172 0 1.4 2
## 4 56 1 1 120 236 0 1 178 0 0.8 2
## 5 57 0 0 120 354 0 1 163 1 0.6 2
## 6 57 1 0 140 192 0 1 148 0 0.4 1
## # ... with 3 more variables: caa <dbl>, thall <dbl>, output <dbl>\[~\]
Analyzing these features, I recognized that the:
- Age -> it refers to the age of the patient
- sex -> it refers if the patient is:
- male (1)
- female (0)
- cp -> it refers to the chest pain type that could be of four types:
- 0 is the typical angina
- 1 is the atypical angina
- 2 is the non-anginal pain
- 3 is the asymptomatic type
- trestbps -> it refers to the resting blood pressure
- chol -> it refers to the serum cholesterol in mg/dl
- fbs -> it refers to the fasting blood sugar that if it is larger than 120 mg/dl is represented with:
- 1 (true)
- 0 (false)
- restecg -> it refers to the resting electrocardiographic results:
- 0 as Normal
- 1 having ST-T wave abnormality (T wave inversions and/or ST elevation or depression of > 0.05 mV)
- showing probable or definite left ventricular hypertrophy by Estes’ criteria
- thalachh -> it refers to the maximum heart rate achieved
- exng -> it refers to the exercise induced angina that could be:
- 1 (yes)
- 0 (no)
- oldpeak -> it refers to the Previous peak
- slp -> it refers to the peak exercise ST segment:
- 0 as up sloping
- 1 as flat
- 2 as down sloping
- caa -> it refers to the number of major vessels that goes from 0 to 3
- thall -> it refers to the maximum heart rate achieved:
- 0 as no-data
- 1 as normal
- 2 as fixed defect
- 3 as reversible defect
- output -> it refers to the fact to have the heart attack or not
\[~\]
The main features detected are:
\[~\]
## age sex cp trtbps
## Min. :29.00 Min. :0.0000 Min. :0.0000 Min. : 94.0
## 1st Qu.:48.00 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:120.0
## Median :55.50 Median :1.0000 Median :1.0000 Median :130.0
## Mean :54.42 Mean :0.6821 Mean :0.9636 Mean :131.6
## 3rd Qu.:61.00 3rd Qu.:1.0000 3rd Qu.:2.0000 3rd Qu.:140.0
## Max. :77.00 Max. :1.0000 Max. :3.0000 Max. :200.0
## chol fbs restecg thalachh
## Min. :126.0 Min. :0.000 Min. :0.0000 Min. : 71.0
## 1st Qu.:211.0 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.:133.2
## Median :240.5 Median :0.000 Median :1.0000 Median :152.5
## Mean :246.5 Mean :0.149 Mean :0.5265 Mean :149.6
## 3rd Qu.:274.8 3rd Qu.:0.000 3rd Qu.:1.0000 3rd Qu.:166.0
## Max. :564.0 Max. :1.000 Max. :2.0000 Max. :202.0
## exng oldpeak slp caa
## Min. :0.0000 Min. :0.000 Min. :0.000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:1.000 1st Qu.:0.0000
## Median :0.0000 Median :0.800 Median :1.000 Median :0.0000
## Mean :0.3278 Mean :1.043 Mean :1.397 Mean :0.7185
## 3rd Qu.:1.0000 3rd Qu.:1.600 3rd Qu.:2.000 3rd Qu.:1.0000
## Max. :1.0000 Max. :6.200 Max. :2.000 Max. :4.0000
## thall output
## Min. :0.000 Min. :0.000
## 1st Qu.:2.000 1st Qu.:0.000
## Median :2.000 Median :1.000
## Mean :2.315 Mean :0.543
## 3rd Qu.:3.000 3rd Qu.:1.000
## Max. :3.000 Max. :1.000\[~\]
The Data Visualization
\[~\]
In this subsection we want to describe graphically the dataset that we are analyzing to have a first taste of what is the better model to use in this case, we want to show below some interesting plots
\[~\]
The visualization using PCA
\[~\]
PCA is a particular tool that allows us to show the behaviour of the patients (treated as points in this space) of our data. We see in this way which are the patients more similar to each other
\[~\]
\[~\]
Its normal to denote that there could be some losses on this representation due to the amount of features that we want to consider.
\[~\]
The Categorical Data
\[~\]
Here, we want to analyze the percentages of each categorical data:
\[~\]
\[~\]
The Quantitative Data
\[~\]
Here, we illustrate the main features of the quantitative data, after have seen the categorical data in the previous section:
\[~\]
hist.and.summary('age', 'Persons Age')## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 29.00 48.00 55.50 54.42 61.00 77.00hist.and.summary('chol', 'Cholestoral')## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 126.0 211.0 240.5 246.5 274.8 564.0hist.and.summary('thalachh', 'Maximum Heart Rate Achieved')## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 71.0 133.2 152.5 149.6 166.0 202.0hist.and.summary('trtbps', 'Resting Blood Pressure')## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 94.0 120.0 130.0 131.6 140.0 200.0hist.and.summary('oldpeak', 'Previous Peak')## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.800 1.043 1.600 6.200\[~\]
After this, we consider also the densities considering the case of having the heart attack and not:
\[~\]
dense.chart('age', 'Persons Age')dense.chart('chol', 'Cholestoral')dense.chart('thalachh', 'Maximum Heart Rate Achieved')dense.chart('trtbps', 'Resting Blood Pressure')dense.chart('oldpeak', 'Previous Peak')The Correlations
\[~\]
Here, we want to highlight the correlations between the features treated in qualitative, quantitative and without any distinction.
\[~\]
The Goal
\[~\]
The main goal is to leverage the main fully Bayesian analysis, based on understanding if a person has a heart attack or not.
So, the response variable \(Y_i\) (that is the “output” feature in our dataset) is the heart attack \(\in \{0,1\}\) and the predictor variables (\(x_i \in \mathbb{R}^{+}\)) have been chosen using the glm() function, used to fit generalized linear models, as we can see below:
\[~\]
##
## Call:
## glm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 +
## x10 + x11 + x12 + x13, family = binomial(), data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5751 -0.3868 0.1514 0.5841 2.6239
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.304169 2.577891 1.282 0.199936
## x1 -0.001469 0.023464 -0.063 0.950062
## x2 -1.750930 0.468199 -3.740 0.000184 ***
## x3 0.847283 0.185548 4.566 4.96e-06 ***
## x4 -0.020188 0.010386 -1.944 0.051916 .
## x5 -0.004489 0.003806 -1.179 0.238252
## x6 0.073463 0.532452 0.138 0.890263
## x7 0.450607 0.348505 1.293 0.196021
## x8 0.023134 0.010449 2.214 0.026835 *
## x9 -0.981017 0.409806 -2.394 0.016672 *
## x10 -0.523604 0.214467 -2.441 0.014630 *
## x11 0.589074 0.349864 1.684 0.092236 .
## x12 -0.826015 0.201922 -4.091 4.30e-05 ***
## x13 -0.887203 0.290730 -3.052 0.002276 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 416.42 on 301 degrees of freedom
## Residual deviance: 210.35 on 288 degrees of freedom
## AIC: 238.35
##
## Number of Fisher Scoring iterations: 6\[~\]
As we can see above, the glm rejects the hypothesis to consider the variables:
- Age ~ \(x_1\)
- Cholestoral ~ \(x_5\)
- fasting blood sugar ~ \(x_6\)
- resting electrocardiographic results ~ \(x_7\)
… and also the intercept isn’t a good choice to admit it in our model.
then, in summary we have N = 302 observations. Then we decided to consider two models and see the difference on which model at the end is the best model of our analysis.
\[~\]
Preliminary brief definitions
\[~\]
Likelihood \(\pi(y_{obs}|\theta)\): measures the goodness of fit of a statistical model to a sample of data for given values of the unknown parameters. It is formed from the joint probability distribution of the sample, but viewed and used as a function of the parameters only, thus treating the random variables as fixed at the observed values.
Logistic regression model: Logistic regression is the appropriate regression analysis to conduct when the dependent variable is dichotomous (binary). Logistic regression is used to describe data and to explain the relationship between one dependent binary variable and one or more nominal, ordinal, interval or ratio-level independent variables.
logit link function: it uses the Cumulative Distribution Function (CDF) of the logistic distribution. A benefit to using Logit is that the coefficients can be interpreted in the terms of odds ratios.
cloglog function: Is unlike Logit and Probit asymmetric. It is best used when the probability of an event is very small or very large. The complementary log-log approaches 0 infinitely slower than any other link function. Cloglog model is closely related to continuous-time models for the occurrence of events, so it has an important application in the area of survival analysis and hazard modeling.
Prior \(\pi(\theta)\): is the probability distribution that would express one’s beliefs about this quantity before some evidence is taken into account.
the deviance: is a goodness-of-fit statistic for a statistical model; it is often used for statistical hypothesis testing. It is a generalization of the idea of using the sum of squares of residuals (RSS) in ordinary least squares to cases where model-fitting is achieved by maximum likelihood.
\[~\]
The First Model (link function with logit)
\[~\]
In the first model, we decide to consider the following logistic regression model with the logit link function:
\[ Y_i \sim Bern(p_i) \\ logit(p_i) = \beta_{2} \cdot x_{2_i} + \beta_{3} \cdot x_{3_i} + \beta_{4} \cdot x_{4_i} +\beta_{8} \cdot x_{8_i} + \beta_{9} \cdot x_{9_i} + \beta_{10} \cdot x_{10_i} + \beta_{11} \cdot x_{11_i} + \beta_{12} \cdot x_{12_i} + \beta_{13} \cdot x_{13_i} \] The prior beta parameters are chosen approximatly considering the estimated averages saw in the previous subsection regarding the outcomes of glm(). So, the beta prior parameters are distributed following a normal distribution as follows:
\[ \beta_2 \sim N(\mu_2 = -1, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_3 \sim N(\mu_3 = 0.5, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_4 \sim N(\mu_4 = 0, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_8 \sim N(\mu_8 = 0, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_9 \sim N(\mu_9 = -0.5, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_{10} \sim N(\mu_{10} = 0.5, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_{11} \sim N(\mu_{11} = 0.5, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_{12} \sim N(\mu_{12} = -0.5, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \beta_{13} \sim N(\mu_{13} = -0.5, \tau^{2} = 1.0 \cdot 10^{-6}) \\ \]
The main reason is that we prefer chosing the binary classification as in this case the logistic regression to model our outcomes based only on 0 or 1.
\[~\]
Implementing RJags
\[~\]
As we can see, we implemented the model using RJags:
\[~\]
- The model is defined as follows:
\[~\]
model <- function(){
# Likelihood
for (i in 1:N){
y[i] ~ dbern(p[i])
logit(p[i]) <- beta2*x2[i] + beta3*x3[i] + beta4*x4[i] + beta8*x8[i] + beta9*x9[i] + beta10*x10[i] + beta11*x11[i] + beta12*x12[i] + beta13*x13[i] # omit intercept and not useful variables
}
# Defining the prior beta parameters
beta2 ~ dnorm(-1, 1.0E-6)
beta3 ~ dnorm(0.5, 1.0E-6)
beta4 ~ dnorm(0, 1.0E-6)
beta8 ~ dnorm(0, 1.0E-6)
beta9 ~ dnorm(-0.5, 1.0E-6)
beta10 ~ dnorm(0.5, 1.0E-6)
beta11 ~ dnorm(0.5, 1.0E-6)
beta12 ~ dnorm(-0.5, 1.0E-6)
beta13 ~ dnorm(-0.5, 1.0E-6)
}remark: that the standard deviation in RJags is considered as the precision, so the \(\tau^2 = \frac{1}{\sigma^{2}}\)
\[~\]
- Data structuring as requested with highlight which parameters to consider:
# Passing the data for RJags
data.jags <- list("y" = y, "N" = N,
"x2" = x2, "x3" = x3, "x4" = x4, "x8" = x8, "x9" = x9, "x10" = x10, "x11" = x11, "x12" = x12, "x13" = x13)
# Defining parameters of interest to show after running RJags
mod.params <- c("beta2", "beta3", "beta4", "beta8", "beta9", "beta10", "beta11", "beta12", "beta13")\[~\]
- Running the model with the following results:
# Run JAGS
# set.seed(123)
n.chains <- 3
mod.fit <- jags(data = data.jags, # DATA
model.file = model, # MODEL
parameters.to.save = mod.params, # TRACKING
n.chains = n.chains, n.iter = 10000, n.burnin = 1000, n.thin=10) # MCMC## module glm loaded## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 302
## Unobserved stochastic nodes: 9
## Total graph size: 3840
##
## Initializing modelmod.fit## Inference for Bugs model at "C:/Users/Jeremy/AppData/Local/Temp/RtmpEhr9h8/model32dc3e3c288f.txt", fit using jags,
## 3 chains, each with 10000 iterations (first 1000 discarded), n.thin = 10
## n.sims = 2700 iterations saved
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta10 -0.488 0.210 -0.913 -0.630 -0.481 -0.342 -0.088 1.001 2700
## beta11 0.685 0.344 -0.004 0.457 0.684 0.917 1.345 1.001 2700
## beta12 -0.830 0.198 -1.238 -0.962 -0.825 -0.698 -0.443 1.000 2700
## beta13 -0.847 0.280 -1.409 -1.032 -0.843 -0.658 -0.311 1.001 2200
## beta2 -1.642 0.445 -2.491 -1.940 -1.633 -1.339 -0.785 1.003 760
## beta3 0.870 0.186 0.517 0.746 0.866 0.995 1.239 1.002 1100
## beta4 -0.014 0.008 -0.031 -0.020 -0.014 -0.009 0.001 1.002 1100
## beta8 0.031 0.008 0.015 0.026 0.031 0.036 0.048 1.003 2700
## beta9 -0.890 0.398 -1.704 -1.155 -0.889 -0.623 -0.136 1.002 1700
## deviance 225.298 4.399 218.687 222.136 224.601 227.809 235.585 1.000 2700
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 9.7 and DIC = 235.0
## DIC is an estimate of expected predictive error (lower deviance is better).\[~\]
In this last running of the MCMC we have interesting results:
mu.vect and sd.vect: these values representing the estimated values after simulating the MCMC.
the percentages (quantiles): represents the quantiles of the interval where the parameter is likely in
RHat: The values near 1 suggest convergence
n.eff: is the effective number of samples to achieve the convergence, the stationary region.
\[~\]
The stationary regions
\[~\]
Here, we show the univariate trace-plots of the simulations of each parameter:
\[~\]
The Cumulative Means
\[~\] In this section we want to see the empirical average of \(\mathbf{\hat{I}}_{t}\) increasing the value of \(t \, = \, 1,...,T\). Before starting, we want to write the formula of \(\mathbf{\hat{I}}_{t}\) that is:
\[~\]
\[ \mathbf{I} \approx \mathbf{\hat{I}}_{t} = \frac{1}{T} \sum_{i=1}^{T} h(\theta^{(i)}) \]
\[~\]
is important to write \(\approx\) because this leverages the SLLN theorem and we want also to confirm this assumption! So, let’s see the results:
\[~\]
\[~\]
Some other interesting plots
\[~\]
We are considering other interesting plots, underlining also each chain alone:
\[~\] Now, we want to see the autocorrelations of each parameter:
## beta10 beta11 beta12 beta13 beta2
## Lag 0 1.00000000 1.000000000 1.000000000 1.0000000000 1.000000000
## Lag 10 -0.01942328 -0.008906691 0.021990698 -0.0027215796 -0.007830619
## Lag 50 -0.01002560 -0.003840567 0.010904558 0.0145067470 0.003030051
## Lag 100 0.00572575 0.020067149 0.002909921 -0.0001148413 0.017380004
## Lag 500 0.03226104 0.015538564 0.023792705 0.0007888182 -0.019429087
## beta3 beta4 beta8 beta9 deviance
## Lag 0 1.0000000000 1.00000000 1.00000000 1.000000000 1.000000000
## Lag 10 0.0344456650 0.02435313 0.04698347 0.001485969 0.029232987
## Lag 50 -0.0008254393 0.02721467 0.00957957 0.011307210 -0.003607147
## Lag 100 -0.0060114153 -0.03236425 0.01009956 -0.039881985 -0.022902722
## Lag 500 0.0243103759 0.03426043 0.03118317 0.011573931 -0.010553928\[~\] The ACFs follow good behaviours, because we need to have samples independent each other during new iterations. Here, we see that the correlations are going further from the first lag, this is strictly decreasing. So, this is a good point!
We can conferm our intuition also with the values returned by autocorr.diag() that it returned the vector of the average autocorrelations across all chains.
\[~\]
The Approximation Errors
\[~\] Here, we want to analyze the best approximation error for the MCMC sampling. We consider, essentially the square root of the MCSE.
The variance formula in the MCMC samplig is:
\[ \mathbf{V}[\hat{I}_{t}] = \frac{Var_{\pi}[h(X_{1})]}{T_{eff}} = \Big( 1 + 2 \sum_{k=1}^{\infty} \rho_{k}\Big)\frac{\sigma^{2}}{T} \] in this case we want to consider the MCSE that is the square root of the formula written above, So the results are (for each chain):
\[~\]
n <- length(colnames(mod.fit$BUGSoutput$sims.matrix))
mcse_dataframe <- data.frame(MCSE_Chain1 = rep(NA, n), MCSE_Chain2 = rep(NA, n), MCSE_Chain3 = rep(NA, n), MCSE_mean = rep(NA, n))
rownames(mcse_dataframe) <- colnames(mod.fit$BUGSoutput$sims.matrix)[1:n]
for(colname in colnames(mod.fit$BUGSoutput$sims.matrix)[1:n]){
for(j in 1:n.chains)
mcse_dataframe[colname, j] <- LaplacesDemon::MCSE(mod.fit$BUGSoutput$sims.array[, j, colname])
mcse_dataframe[colname, "MCSE_mean"] <- mean(unlist(mcse_dataframe[colname, c(1:3)]))
}
mcse_dataframe## MCSE_Chain1 MCSE_Chain2 MCSE_Chain3 MCSE_mean
## beta10 0.0066750908 0.0070902007 0.0067937566 0.0068530160
## beta11 0.0110584487 0.0122883577 0.0116360858 0.0116609641
## beta12 0.0065767369 0.0078695947 0.0067823288 0.0070762201
## beta13 0.0093975655 0.0090483166 0.0094719588 0.0093059470
## beta2 0.0146390539 0.0157325177 0.0145421558 0.0149712425
## beta3 0.0069944633 0.0060050009 0.0066680984 0.0065558542
## beta4 0.0003077438 0.0002888290 0.0002913331 0.0002959686
## beta8 0.0002753919 0.0003032742 0.0003012428 0.0002933030
## beta9 0.0130918777 0.0133971276 0.0135249396 0.0133379817
## deviance 0.1506761221 0.1738054127 0.1517291188 0.1587368845\[~\]
Hypothesis and Testing
\[~\]
\[~\]
The Predictions
\[~\]
The Correlations
\[~\] Here, we bring the correlations between all the values calculated during the MCMC sampling, considering the joining matrix:
\[~\]
The Second Model
\[~\]
The Cumulative Means
\[~\]
The Approximation Errors
\[~\]
The Predictions
\[~\]
The Correlations
\[~\]
The comparison between these two models
\[~\]
The Conclusions
\[~\]
\[~\]
The References
\[~\]